########################################################################################################
##                                                                                                    ##
##  This is the replication file for the real data analysis in Tsai, Tsung-han. 2014.                 ##
##  "A Bayesian Approach to Dynamic Panel Data Models with Rarely Changing Variables."                ##
##                                                                                                    ##
########################################################################################################
## Before running, change the working directory in line 26 to the appropriate one.
##
## Require the following files:
## (1) panel_functions.R: auxiliary functions
## (2) bayes.mlm.R: The syntax for the Bayesian multilevel model
## (3) bayes.endo4.R: The syntax for the Bayesian simultaneous equation model (BSEM)
## (4) LAimplag1.csv: Latin American data
##
## Require the following packages:
## (1) plm
## (2) arm
## (3) R2jags
## (4) pcse
## (5) MCMCpack

## CLEAN UP
rm(list=ls());

## SET WORKING DIRECTORY
setwd("");

## LOAD REQUIRED FUNCTIONS
source("panel_functions.R");

## LOAD PACKAGES.
library(ecodist);  # GENERATE CORRELATED SAMPLES: corgen
library(plm);
library(arm);
library(R2jags);
library(pcse);
library(MCMCpack);

## SET THE SEED, WHICH TELLS US WHERE THE RANDOM PROCESS STARTS.
set.seed(04102014);


###################################################################################
## Data Analysis
###################################################################################
## Load the data set
IMP.Huberlag1 <- read.csv("LAimplag1.csv", header=T);
#dim(IMP.Huberlag1);

################################################################
## MODEL 1: POOLED OLS WITH PCSE
pooled.ols <- lm(log.csssw ~ l1.log.csssw + dempolity + rgdpch.1000 + ubnpop.per + pop65wdi + trade + fdigdp, data=IMP.Huberlag1);
round(summary(pooled.ols)$coef, 3);
pols.pcse <- pcse(pooled.ols, groupN=IMP.Huberlag1$country, groupT=IMP.Huberlag1$year);
summary(pols.pcse);

coef.ols <- round(cbind(pols.pcse$b, pols.pcse$pcse, pols.pcse$b-1.64*pols.pcse$pcse, pols.pcse$b+ 1.64*pols.pcse$pcse), 3);
colnames(coef.ols) <- c("Est", "se", "5%", "95%");

################################################################
## MODEL 2: FIXED EFFECTS MODEL
fe.mod <- plm(log.csssw ~ l1.log.csssw + dempolity + rgdpch.1000 + ubnpop.per + pop65wdi + trade + fdigdp, data=IMP.Huberlag1, index=c("country", "year"), model="within"); 
round(summary(fe.mod)$coef, 3);

coef.fe <- round(cbind(summary(fe.mod)$coef[,1], summary(fe.mod)$coef[,2], summary(fe.mod)$coef[,1]-1.96*summary(fe.mod)$coef[,2], summary(fe.mod)$coef[,1]+1.96*summary(fe.mod)$coef[,2]), 3);
coef.fe <- rbind(0, coef.fe);
colnames(coef.fe) <- c("Est", "se", "5%", "95%");

#write.table(coef.fe, file="coef_fe.txt", sep=",", row.names=TRUE, col.names=TRUE);
coef.fe <- read.table(file="coef_fe.txt", header=TRUE, sep=",");

################################################################
## MODEL 3: RANDOM EFFECTS MODEL
re.mod <- plm(log.csssw ~ l1.log.csssw + dempolity + rgdpch.1000 + ubnpop.per + pop65wdi + trade + fdigdp, data=IMP.Huberlag1, index=c("country", "year"), model="random", random.method="amemiya"); 
round(summary(re.mod)$coef, 3);

coef.re <- round(cbind(summary(re.mod)$coef[,1], summary(re.mod)$coef[,2], summary(re.mod)$coef[,1]-1.96*summary(re.mod)$coef[,2], summary(re.mod)$coef[,1]+1.96*summary(re.mod)$coef[,2]), 3);
colnames(coef.re) <- c("Est", "se", "2.5%", "97.5%");

#write.table(coef.re, file="coef_re.txt", sep=",", row.names=TRUE, col.names=TRUE);
coef.re <- read.table(file="coef_re.txt", header=TRUE, sep=",");

################################################################
## MODEL 4: FIXED EFFECTS VECTOR DECOMPOSITION (fevd) MODEL
# FIRST STAGE: FIT A STANDARD FIXED EFFECTS MODEL
fevd.1 <- plm(log.csssw ~ l1.log.csssw + dempolity + rgdpch.1000 + ubnpop.per + pop65wdi + trade + fdigdp, data=IMP.Huberlag1, index=c("country", "year"), model="within"); 
#summary(fevd.1)$coef[,1:2];

b <- coef(fevd.1); b <- matrix(b, nrow=length(b), ncol=1, byrow=FALSE);
X.bar <- IMP.Huberlag1[, c("ctry.l1logcsssw","ctry.polity","ctry.rgdp1000","ctry.ubnpop.per","ctry.pop65wdi","ctry.trade","ctry.fdigdp")];
X.bar <- as.matrix(X.bar);
e <- residuals(fevd.1);

e <- matrix(e, nrow=18, ncol=30, byrow=TRUE);
ebar <- rep(rowMeans(e), each=30);
ebar <- matrix(ebar, nrow=18*30, ncol=1, byrow=FALSE);

u.hat <- IMP.Huberlag1$ctry.logcsssw - X.bar %*% b - ebar;

# SECOND STAGE: REGRESS THE ESTIMATED UNIT EFFECTS ON THE RARELY CHANGING VARIABLES
fevd.2 <- lm(u.hat ~ -1 + IMP.Huberlag1$dempolity);
#summary(fevd.2)$coef[,c(1,2)];

IMP.Huberlag1$h <- residuals(fevd.2);

# THIRD STAGE: RERUN THE FULL MODEL WITHOUT UNIT EFFECTS BUT INCLUDE h
fevd.3 <- lm(log.csssw ~ l1.log.csssw + dempolity + rgdpch.1000 + ubnpop.per + pop65wdi + trade + fdigdp+h, data=IMP.Huberlag1);
#summary(fevd.3)$coef[2:8,1];

round(summary(fevd.3)$coef, 3);

coef.fevd <- round(cbind(summary(fevd.3)$coef[,1], summary(fevd.3)$coef[,2], summary(fevd.3)$coef[,1]-1.96*summary(fevd.3)$coef[,2], summary(fevd.3)$coef[,1]+1.96*summary(fevd.3)$coef[,2]), 3);
colnames(coef.fevd) <- c("Est", "se", "2.5%", "97.5%");

#write.table(coef.fevd, file="coef_fevd.txt", sep=",", row.names=TRUE, col.names=TRUE);
coef.fevd <- read.table(file="coef_fevd.txt", header=TRUE, sep=",");

################################################################
## MODEL 5: ML MULTILEVEL MODEL WITHOUT GROUP MEANS
ADL11.mlm1 <- lmer(log.csssw ~ l1.log.csssw + dempolity + rgdpch.1000 + ubnpop.per + pop65wdi + trade + fdigdp + (1 | country), data=IMP.Huberlag1);

summary(ADL11.mlm1);

################################################################
## MODEL 6: ML MULTILEVEL MODEL WITH GROUP MEANS
ADL11.mlm2 <- lmer(log.csssw ~ l1.log.csssw + dempolity + rgdpch.1000 + ubnpop.per + pop65wdi + trade + fdigdp + ctry.polity + ctry.rgdp1000 + ctry.ubnpop.per + ctry.pop65wdi + ctry.trade + ctry.fdigdp + (1 | country), data=IMP.Huberlag1);

summary(ADL11.mlm2)$coef;

coef.mlm <- round(cbind(summary(ADL11.mlm2)$coef[1:8,1:2], summary(ADL11.mlm2)$coef[1:8,1]-1.64*summary(ADL11.mlm2)$coef[1:8,2], summary(ADL11.mlm2)$coef[1:8,1]+ 1.64*summary(ADL11.mlm2)$coef[1:8,2]), 3);
colnames(coef.mlm) <- c("Est", "se", "5%", "95%");

################################################################
## MODEL 7: BAYESIAN MULTILEVEL MODEL WITH GROUP MEANS
N <- dim(IMP.Huberlag1)[1];
J <- length(unique(IMP.Huberlag1$country));
country <- as.numeric(IMP.Huberlag1$country);

y <- IMP.Huberlag1$log.csssw;

attach(IMP.Huberlag1);

# COVARIATES OF Y REGRESSION
X.1 <- cbind(l1.log.csssw, dempolity, rgdpch.1000, ubnpop.per, pop65wdi, trade, fdigdp);
M.1 <- dim(X.1)[2];

# COVARIATES OF THE SECOND LEVEL REGRESSION
Z.1 <- round(cbind(1, unique(ctry.polity), unique(ctry.rgdp1000), unique(ctry.pop65wdi), unique(ctry.trade), unique(ctry.fdigdp)), 3);
L <- dim(Z.1)[2];

detach(IMP.Huberlag1);

data1 <- list("y", "N", "J", "country", "X.1", "M.1", "Z.1", "L");

# INITIAL VALUES
inits1 <- list(
	list( tau.y=1, B=rep(0, M.1), tau.d=0.1, G=rep(-3, L) ),
	list( tau.y=0.1, B=rep(3, M.1), tau.d=0.5, G=rep(3, L) ),
	list( tau.y=0.5, B=rep(-3, M.1), tau.d=1, G=rep(0, L) )
	)

parameters1 <- c("B", "d", "G");

n.chains <- length(inits1)
n.iteration <- 50000;
frac.discard <- 0.5;  # FRACTION OF SAMPLES DISCARDED
n.discard <- n.iteration*frac.discard; n.discard;  # BURNIN-IN PERIOD
n.thining <- 10;
n.saved.chin <- (n.iteration - n.discard)/n.thining; n.saved.chin;  # SAMPLES AFTER THINING SAVED FOR EACH CHAIN
n.saved.chin*n.chains;  # TOTAL SAMPLES SAVED


## hand-off to JAGS
bmlm <- jags(model.file="bayes.mlm.R",
             data = data1,
             inits = inits1,
             parameters.to.save = parameters1,
             n.chains = n.chains,   # (DEFAULT: 1)
             n.iter = n.iteration,
             n.burnin = n.discard,
             n.thin = n.thining,
             DIC = TRUE
             );

print(bmlm);

##
mcmcburnin.mcmclist <- as.mcmc(bmlm);


## REPORT THE RESULTS
out.bmlm <- as.matrix(mcmcburnin.mcmclist);
dim(out.bmlm);

est.bmlm <- round(cbind(apply(out.bmlm, 2, mean), apply(out.bmlm, 2, sd), t(apply(out.bmlm, 2, quantile, prob=c(0.05, 0.95))) ), 3);
colnames(est.bmlm)[1:2] <- c("mean", "sd");
coef.bmlm <- est.bmlm[c(27, 1:7),];

       mean    sd     5%    95%
G[1] -1.867 0.565 -2.807 -0.988
B[1]  0.450 0.100  0.389  0.512
B[2]  0.013 0.030  0.002  0.024
B[3] -0.076 0.096 -0.142 -0.014
B[4]  0.016 0.033  0.003  0.028
B[5]  0.039 0.159 -0.076  0.154
B[6] -0.002 0.021 -0.005  0.001
B[7]  0.001 0.037 -0.025  0.026


################################################################
## MODEL 8: BAYESIAN SIMULTANEOUS EQUATION MODEL
## SET THE DATA
N <- dim(IMP.Huberlag1)[1];
J <- length(unique(IMP.Huberlag1$country));
country <- IMP.Huberlag1$country;

y <- IMP.Huberlag1$log.csssw;

attach(IMP.Huberlag1);

# COVARIATES OF Y REGRESSION
X.1 <- cbind(1, l1.log.csssw, dempolity, rgdpch.1000, ubnpop.per, pop65wdi, trade, fdigdp);
M.1 <- dim(X.1)[2];

w1 <- dempolity;
w2 <- rgdpch.1000;
w3 <- ubnpop.per;
w4 <- pop65wdi;
w5 <- trade;
w6 <- fdigdp;

# COVARIATES OF w1 REGRESSION: dempolity
W.1 <- matrix(rep(1,N), ncol=1);
L.1 <- dim(W.1)[2];

# COVARIATES OF w2 REGRESSION: rgdpch.1000
W.2 <- matrix(rep(1,N), ncol=1);
L.2 <- dim(W.2)[2];

# COVARIATES OF w3 REGRESSION: ubnpop.per
W.3 <- matrix(rep(1,N), ncol=1);
L.3 <- dim(W.3)[2];

# COVARIATES OF w4 REGRESSION: pop65wdi
W.4 <- matrix(rep(1,N), ncol=1);
L.4 <- dim(W.4)[2];

# COVARIATES OF w5 REGRESSION: trade
W.5 <- matrix(rep(1,N), ncol=1);
L.5 <- dim(W.5)[2];

# COVARIATES OF w6 REGRESSION: fdigdp
W.6 <- matrix(rep(1,N), ncol=1);
L.6 <- dim(W.6)[2];

detach(IMP.Huberlag1);

K <- 7; # THE NUMBER OF ENDOGENOUS REGRESSIONS PLUS UNIT EFFECTS

W <- diag(K);

data1 <- list("y", "N", "J", "country", "X.1", "M.1", "K", "W", "w1", "W.1", "L.1", "w2", "W.2", "L.2", "w3", "W.3", "L.3", "w4", "W.4", "L.4", "w5", "W.5", "L.5", "w6", "W.6", "L.6");

# INITIAL VALUES FROM MLE
inits1 <- list(
	list(tau.y=1, B=rep(0, M.1), D=array(rep(0,J*K), c(J,K)), Tau.D=W, G.1=rep(0, L.1), G.2=rep(0, L.2), G.3=rep(0, L.3), G.4=rep(0, L.4), G.5=rep(0, L.5), G.6=rep(0, L.6), tau.w1=1, tau.w2=1, tau.w3=1, tau.w4=1, tau.w5=1, tau.w6=1 ),
	list(tau.y=.1, B=rep(3, M.1), D=array(rep(3,J*K), c(J,K)), Tau.D=W, G.1=rep(3, L.1), G.2=rep(-3, L.2), G.3=rep(3, L.3), G.4=rep(-3, L.4), G.5=rep(-3, L.5), G.6=rep(3, L.6), tau.w1=.1, tau.w2=.1, tau.w3=.1, tau.w4=.1, tau.w5=.1, tau.w6=.1 ),
	list(tau.y=.5, B=rep(-3, M.1), D=array(rep(-3,J*K), c(J,K)), Tau.D=W, G.1=rep(-3, L.1), G.2=rep(3, L.2), G.3=rep(-3, L.3), G.4=rep(3, L.4), G.5=rep(3, L.5), G.6=rep(-3, L.6), tau.w1=.5, tau.w2=.5, tau.w3=.5, tau.w4=.5, tau.w5=.5, tau.w6=.5 )
	);


parameters1 <- c("sigma.y", "B", "sigma.d", "sigma.w1", "sigma.w2", "sigma.w3", "sigma.w4", "sigma.w5", "sigma.w6", "sigma.22", "sigma.33", "sigma.44", "sigma.55", "sigma.66", "sigma.77", "rho.12", "rho.13", "rho.14", "rho.15", "rho.16", "rho.17");

n.chains <- length(inits1);
n.iteration <- 50000;
frac.discard <- 0.5;  # FRACTION OF SAMPLES DISCARDED
n.discard <- n.iteration*frac.discard; n.discard;  # BURNIN-IN PERIOD
n.thining <- 5;
n.saved.chin <- (n.iteration - n.discard)/n.thining; n.saved.chin;  # SAMPLES AFTER THINING SAVED FOR EACH CHAIN
n.saved.chin*n.chains;  # TOTAL SAMPLES SAVED


## hand-off to JAGS
bml1 <- jags(model.file="bayes.endo4.R",
             data = data1,
             inits = inits1,
             parameters.to.save = parameters1,
             n.chains = n.chains,   # (DEFAULT: 1)
             n.iter = n.iteration,
             n.burnin = n.discard,
             n.thin = n.thining,
             DIC = TRUE
             );

print(bml1);


##
mcmcburnin.mcmclist <- as.mcmc(bml1);


## REPORT THE RESULTS
out <- as.matrix(mcmcburnin.mcmclist);
dim(out);

write.csv(out, "LA_bayes.csv", row.names=FALSE);



##############################################################################
## LOAD THE MCMC OUTPUT
out <- read.csv("LA_bayes.csv");

est.bayes <- round(cbind(apply(out, 2, mean), apply(out, 2, sd), t(apply(out, 2, quantile, prob=c(0.05, 0.95))) ), 3);
colnames(est.bayes)[1:2] <- c("mean", "sd");
coef.bayes <- est.bayes[1:8,];


##  PRESENT A PLOT WITH MEANS AND 95% CREDIBLE INTERVALS
coef.ols;
coef.fe;
coef.fevd;
coef.re;
coef.mlm;
coef.bayes;


mod.label <- c("OLS", "FE", "FEVD", "RE", "MLM", "BSEM");
loc <- seq(1.5,5.5,length.out=6);

par(mfrow=c(1,1), mar=c(4,4,2,0), oma=c(1,2,1,1));

#######################  Estimates of Intercept  #####################
intercept <- rbind(coef.ols[1,], coef.fe[1,], coef.fevd[1,], coef.re[1,], coef.mlm[1,], coef.bayes[1,]);
rownames(intercept) <- c("ols", "fe", "fevd", "re", "mlm", "bsem");

plot(1:6 ~ intercept[1:6,1], xlim=c(min(intercept), max(intercept)), type="n", xlab="(a) Estimates of the Intercept", ylab="", main="", axes=FALSE);

# X-axis
axis(side=1, las=1, tick=TRUE, cex.axis=1);
axis(side=3, las=1, tick=TRUE, cex.axis=1);

# Y-axis
axis(side=2, at=loc, label=mod.label, las=1, tick=FALSE, , cex.axis=0.8);
abline(h=loc, lty=2, lwd=0.5, col="light gray");

# 90% intervals
segments(intercept[,3], loc, intercept[,4], loc, col="grey50", lwd=3);

# Means
points(x=intercept[,1], y=loc, pch="+", cex=1.5);


#######################  Estimates of LDV  #####################
ldv <- rbind(coef.ols[2,], coef.fe[2,], coef.fevd[2,], coef.re[2,], coef.mlm[2,], coef.bayes[2,]);
rownames(ldv) <- c("ols", "fe", "fevd", "re", "mlm", "bsem");

plot(1:6 ~ ldv[1:6,1], xlim=c(0.3, 0.8), type="n", xlab="Estimates of the LDV", ylab="", main="", axes=FALSE, cex.lab=1.2);

# X-axis
axis(side=1, las=1, tick=TRUE, cex.axis=1.4);
axis(side=3, las=1, tick=TRUE, cex.axis=1.4);

# Y-axis
axis(side=2, at=loc, label=mod.label, las=1, tick=FALSE, cex.axis=1.3);
abline(h=loc, lty=2, lwd=0.5, col="light gray");

# Means
points(x= ldv[,1], y=loc, pch="+", cex=2.5);

# 90% intervals
segments(ldv[,3], loc, ldv[,4], loc, col="grey50", lwd=4);


#######################  Estimates of Democracy  #####################
dem <- rbind(coef.ols[3,], coef.fe[3,], coef.fevd[3,], coef.re[3,], coef.mlm[3,], coef.bayes[3,]);
rownames(dem) <- c("ols", "fe", "fevd", "re", "mlm", "bsem");

plot(1:6 ~ dem[1:6,1], xlim=c(-0.04, 0.04), type="n", xlab="(a) Estimates of Democracy", ylab="", main="", axes=FALSE, cex.lab=1.2);

# X-axis
axis(side=1, las=1, tick=TRUE, cex.axis=1.4);
axis(side=3, las=1, tick=TRUE, cex.axis=1.4);

# Y-axis
axis(side=2, at=loc, label=mod.label, las=1, tick=FALSE, cex.axis=1.3);
abline(h=loc, lty=2, lwd=0.5, col="light gray");

#
abline(v=0, lty=1, lwd=1, col="grey80");

# Means
points(x= dem[,1], y=loc, pch="+", cex=2.5);

# 90% intervals
segments(dem[,3], loc, dem[,4], loc, col="grey50", lwd=4);


#######################  Estimates of GDPPC  #####################
gdppc <- rbind(coef.ols[4,], coef.fe[4,], coef.fevd[4,], coef.re[4,], coef.mlm[4,], coef.bayes[4,]);
rownames(gdppc) <- c("ols", "fe", "fevd", "re", "mlm", "bsem");

plot(1:6 ~ gdppc[1:6,1], xlim=c(-0.2, 0.1), type="n", xlab="(b) Estimates of GDPPC", ylab="", main="", axes=FALSE, cex.lab=1.2);

# X-axis
axis(side=1, las=1, tick=TRUE, cex.axis=1.4);
axis(side=3, las=1, tick=TRUE, cex.axis=1.4);

# Y-axis
axis(side=2, at=loc, label=mod.label, las=1, tick=FALSE, cex.axis=1.3);
abline(h=loc, lty=2, lwd=0.5, col="light gray");

#
abline(v=0, lty=1, lwd=1, col="grey80");

# Means
points(x= gdppc[,1], y=loc, pch="+", cex=2.5);

# 90% intervals
segments(gdppc[,3], loc, gdppc[,4], loc, col="grey50", lwd=4);


#######################  Estimates of Urban Population  #####################
ubnpop <- rbind(coef.ols[5,], coef.fe[5,], coef.fevd[5,], coef.re[5,], coef.mlm[5,], coef.bayes[5,]);
rownames(ubnpop) <- c("ols", "fe", "fevd", "re", "mlm", "bsem");

plot(1:6 ~ ubnpop[1:6,1], xlim=c(-0.01, 0.04), type="n", xlab="(c) Estimates of Urban Population", ylab="", main="", axes=FALSE, cex.lab=1.2);

# X-axis
axis(side=1, las=1, tick=TRUE, cex.axis=1.4);
axis(side=3, las=1, tick=TRUE, cex.axis=1.4);

# Y-axis
axis(side=2, at=loc, label=mod.label, las=1, tick=FALSE, cex.axis=1.3);
abline(h=loc, lty=2, lwd=0.5, col="light gray");

#
abline(v=0, lty=1, lwd=1, col="grey80");

# Means
points(x= ubnpop[,1], y=loc, pch="+", cex=2.5);

# 90% intervals
segments(ubnpop[,3], loc, ubnpop[,4], loc, col="grey50", lwd=4);


#######################  Estimates of Aged Population  #####################
pop65 <- rbind(coef.ols[6,], coef.fe[6,], coef.fevd[6,], coef.re[6,], coef.mlm[6,], coef.bayes[6,]);
rownames(pop65) <- c("ols", "fe", "fevd", "re", "mlm", "bsem");

plot(1:6 ~ pop65[1:6,1], xlim=c(-0.2, 0.2), type="n", xlab="(d) Estimates of Aged Population", ylab="", main="", axes=FALSE, cex.lab=1.2);

# X-axis
axis(side=1, las=1, tick=TRUE, cex.axis=1.4);
axis(side=3, las=1, tick=TRUE, cex.axis=1.4);

# Y-axis
axis(side=2, at=loc, label=mod.label, las=1, tick=FALSE, cex.axis=0.8);
abline(h=loc, lty=2, lwd=0.5, col="light gray");

#
abline(v=0, lty=1, lwd=1, col="grey80");

# Means
points(x= pop65[,1], y=loc, pch="+", cex=2.5);

# 90% intervals
segments(pop65[,3], loc, pop65[,4], loc, col="grey50", lwd=4);


#######################  Estimates of Trade  #####################
trade <- rbind(coef.ols[7,], coef.fe[7,], coef.fevd[7,], coef.re[7,], coef.mlm[7,], coef.bayes[7,]);
rownames(trade) <- c("ols", "fe", "fevd", "re", "mlm", "bsem");

plot(1:6 ~ trade[1:6,1], xlim=c(-0.01, 0.01), type="n", xlab="(e) Estimates of Trade", ylab="", main="", axes=FALSE, cex.lab=1.2);

# X-axis
axis(side=1, las=1, tick=TRUE, cex.axis=1.4);
axis(side=3, las=1, tick=TRUE, cex.axis=1.4);

# Y-axis
axis(side=2, at=loc, label=mod.label, las=1, tick=FALSE, cex.axis=1.3);
abline(h=loc, lty=2, lwd=0.5, col="light gray");

#
abline(v=0, lty=1, lwd=1, col="grey80");

# Means
points(x= trade[,1], y=loc, pch="+", cex=2.5);

# 90% intervals
segments(trade[,3], loc, trade[,4], loc, col="grey50", lwd=4);


#######################  Estimates of FDI  #####################
fdi <- rbind(coef.ols[8,], coef.fe[8,], coef.fevd[8,], coef.re[8,], coef.mlm[8,], coef.bayes[8,]);
rownames(fdi) <- c("ols", "fe", "fevd", "re", "mlm", "bsem");

plot(1:6 ~ trade[1:6,1], xlim=c(-0.1, 0.1), type="n", xlab="(f) Estimates of FDI", ylab="", main="", axes=FALSE, cex.lab=1.2);

# X-axis
axis(side=1, las=1, tick=TRUE, cex.axis=1.4);
axis(side=3, las=1, tick=TRUE, cex.axis=1.4);

# Y-axis
axis(side=2, at=loc, label=mod.label, las=1, tick=FALSE, cex.axis=1.3);
abline(h=loc, lty=2, lwd=0.5, col="light gray");

#
abline(v=0, lty=1, lwd=1, col="grey80");

# Means
points(x= fdi[,1], y=loc, pch="+", cex=2.5);

# 90% intervals
segments(fdi[,3], loc, fdi[,4], loc, col="grey50", lwd=4);



############################################################################################
# SINCE THE DISTRIBUTIONS CORRELATION ESTIMATES ARE NOT SYMMETRIC, HIGHEST POSTERIOR DENSITY (HPD) INTERVALS ARE CALCULATED
## USE "HPDint" FUNCTION
w1 <- dempolity;
w2 <- rgdpch.1000;
w3 <- ubnpop.per;
w4 <- pop65wdi;
w5 <- trade;
w6 <- fdigdp;


corr.spend <- matrix(NA, nrow=6, ncol=6);
colnames(corr.spend) <- c("mean", "sd", "90% HPD Lower", "90% HPD Upper", "80% HPD Lower", "80% HPD Upper");
rownames(corr.spend) <- c("rho12", "rho13", "rho14", "rho15", "rho16", "rho17");

corr.spend[1,1:2] <- est.bayes["rho.12",1:2];
corr.spend[1,3:4] <- HPDint(out[,"rho.12"], prob=0.90)$HPDinterval;
corr.spend[1,5:6] <- HPDint(out[,"rho.12"], prob=0.80)$HPDinterval;

corr.spend[2,1:2] <- est.bayes["rho.13",1:2];
corr.spend[2,3:4] <- HPDint(out[,"rho.13"], prob=0.90)$HPDinterval;
corr.spend[2,5:6] <- HPDint(out[,"rho.13"], prob=0.80)$HPDinterval;

corr.spend[3,1:2] <- est.bayes["rho.14",1:2];
corr.spend[3,3:4] <- HPDint(out[,"rho.14"], prob=0.90)$HPDinterval;
corr.spend[3,5:6] <- HPDint(out[,"rho.14"], prob=0.80)$HPDinterval;

corr.spend[4,1:2] <- est.bayes["rho.15",1:2];
corr.spend[4,3:4] <- HPDint(out[,"rho.15"], prob=0.90)$HPDinterval;
corr.spend[4,5:6] <- HPDint(out[,"rho.15"], prob=0.80)$HPDinterval;

corr.spend[5,1:2] <- est.bayes["rho.16",1:2];
corr.spend[5,3:4] <- HPDint(out[,"rho.16"], prob=0.90)$HPDinterval;
corr.spend[5,5:6] <- HPDint(out[,"rho.16"], prob=0.80)$HPDinterval;

corr.spend[6,1:2] <- est.bayes["rho.17",1:2];
corr.spend[6,3:4] <- HPDint(out[,"rho.17"], prob=0.90)$HPDinterval;
corr.spend[6,5:6] <- HPDint(out[,"rho.17"], prob=0.80)$HPDinterval;

corr.spend <- corr.spend[6:1,];

## PLOT CORRELATIONS

######################################################################################################
##  Figure                                                                                          ##
######################################################################################################
##  PRESENT A PLOT WITH MEANS AND 90% CREDIBLE INTERVALS
var.label <- c("DEM", "GDPPC", "UBNPOP", "POP65", "TRADE", "FDI");

y.pos <- seq(from=1.5, to=5.5, length=6);


par(mfrow=c(1,1), mar=c(2,4,6,2), oma=c(4,3,0,0));

plot(1:6 ~ corr.spend[,1], type="n", xlab="", ylab="", main="", axes=FALSE, xlim=c(-1,1));

# X-axis
axis(side=1, las=1, tick=TRUE, cex.axis=1.3);
axis(side=3, las=1, tick=TRUE, cex.axis=1.3);

# Y-axis
axis(side=2, at= y.pos, label=rev(var.label), las=1, tick=FALSE, , cex.axis=1.2);
abline(h= y.pos, lty=2, lwd=0.5, col="light gray");

# 90% credible intervals
segments(corr.spend[,3], y.pos, corr.spend[,4], y.pos, col="grey50", lwd=3);

# 80% credible intervals
segments(corr.spend[,5], y.pos, corr.spend[,6], y.pos, col="grey10", lwd=5)

# Means
points(x= corr.spend[,1], y= y.pos, pch="+", cex=2);



#################################################
## LONG-TERM EFFECTS OF REGIMES
long.run <- out[,3]/(1-out[,2]);
quantile(long.run, prob=c(0.05,0.95));
mean(long.run);
sd(long.run);
HPDint(long.run, prob=0.8);

hist(long.run, breaks=30, probability=TRUE, xlim=c(-.5,.5), main="", xlab="Long-term Effect", ylab="", col="grey90", border="black", axes=FALSE, cex.lab=1.5);
axis(side=1, at=seq(-1, 1, .2), tick=TRUE, cex.axis=1.5);
rug(c(-0.1, 0.181), col="red", lwd=5);



## STEP RESPONSE FUNCTION
time0 <- quantile(out[,3], probs=c(.05, .5, .95))
time1 <- quantile(out[,3] + out[,3]*out[,2], probs=c(.05, .5, .95))
time2 <- quantile(out[,3] + out[,3]*out[,2] + out[,3]*out[,2]^2, probs=c(.05, .5, .95))
time3 <- quantile(out[,3] + out[,3]*out[,2] + out[,3]*out[,2]^2 + out[,3]*out[,2]^3, probs=c(.05, .5, .95))
time4 <- quantile(out[,3] + out[,3]*out[,2] + out[,3]*out[,2]^2 + out[,3]*out[,2]^3 + out[,3]*out[,2]^4, probs=c(.05, .5, .95))
time5 <- quantile(out[,3] + out[,3]*out[,2] + out[,3]*out[,2]^2 + out[,3]*out[,2]^3 + out[,3]*out[,2]^4 + out[,3]*out[,2]^5, probs=c(.05, .5, .95))


o.change <- rbind(time0,time1,time2,time3,time4,time5)
time.per <- c(0:5)


par(mfrow=c(1,1), mar=c(4,4,2,0), oma=c(1,2,1,0));

plot(y=o.change[,2], x=time.per, ylim=c(-.01, .05), type="b", pch="+", ylab="Dynamic Effects", xlab="Time Points", col="black", lwd=2.5, main="", axes=FALSE, cex.lab=1.3)

abline(h=0, col="red")
segments(time.per, o.change[,1], time.per, o.change[,3], col="grey70")
points(y=o.change[,2], x=time.per, pch="+", cex=1.5)
axis(side=1, at=seq(0, 5), tick=TRUE, cex.axis=1.2)
axis(side=2, at=seq(-.01, .05, .01), tick=TRUE, cex.axis=1.2, las=1)



